12 Session 12: Area Data IV

Practice questions

Answer the following questions:

1. How are row-standardized and binary spatial weights interpreted?
2. What is the reason for using a Bonferroni correction for multiple tests?
3. What types of spatial patterns can the local version of Moran’s I detect?
4. What types of spatial patterns can the \(G_i(d)\) statistic detect?
5. What is the utility of detecting hot and cold spatial spots?

Learning objectives

In this activity, you will:

1. Calculate Moran’s I coefficient of autocorrelation for area data.
2. Create Moran’s scatterplots.
3. Examine the results of the tests/scatterplots for further insights.
4. Think about ways to decide whether a landscape is random when working with area data.

Suggested reading

O’Sullivan D and Unwin D (2010) Geographic Information Analysis, 2nd Edition, Chapter 7. John Wiley & Sons: New Jersey.

Preliminaries

For this activity you will need the following:

• This R markdown notebook.
• A shape file called Hamilton CMA CT.

It is good practice to clear the working space to make sure that you do not have extraneous items there when you begin your work. The command in R to clear the workspace is rm (for “remove”), followed by a list of items to be removed. To clear the workspace from all objects, do the following:

rm(list = ls())

Note that ls() lists all objects currently on the worspace.

Load the libraries you will use in this activity. In addition to tidyverse, you will need spdep, a package designed for the analysis of spatial data (you can learn about spdep here and here):

library(tidyverse)

-- Attaching packages --------------------------------- tidyverse 1.2.1 --
v ggplot2 2.2.1     v purrr   0.2.4
v tibble  1.4.2     v dplyr   0.7.4
v tidyr   0.8.0     v stringr 1.2.0
v readr   1.1.1     v forcats 0.2.0
-- Conflicts ------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(rgdal)

Loading required package: sp
rgdal: version: 1.2-16, (SVN revision 701)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
 Path to GDAL shared files: C:/Users/Antonio/Documents/R/win-library/3.4/rgdal/gdal
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: C:/Users/Antonio/Documents/R/win-library/3.4/rgdal/proj
 Linking to sp version: 1.2-7 

library(broom)
library(spdep)

Loading required package: Matrix

Attaching package: <U+393C><U+3E31>Matrix<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:tidyr<U+393C><U+3E32>:

    expand

Loading required package: spData

Begin by loading the data that you will use in this activity:

Hamilton_CT <- readOGR(".", layer = "Hamilton CMA CT", integer64 = "allow.loss")

OGR data source with driver: ESRI Shapefile 
Source: ".", layer: "Hamilton CMA CT"
with 188 features
It has 255 fields
Integer64 fields read as signed 32-bit integers:  ID POPULATION PRIVATE_DW OCCUPIED_D ALL_AGES AGE_4 AGE_5_TO_9 AGE_10_TO_ AGE_15_TO_ AGE_15 AGE_16 AGE_17 AGE_18 AGE_19 AGE_20_TO_ AGE_25_TO_ AGE_30_TO_ AGE_35_TO_ AGE_40_TO_ AGE_45_TO_ AGE_50_TO_ AGE_55_TO_ AGE_60_TO_ AGE_65_TO_ AGE_70_TO_ AGE_75_TO_ AGE_80_TO_ AGE_85 MEDIAN_AGE MALE_ALL_A MALE_4 MALE_5_TO_ MALE_10_TO MALE_15_TO MALE_15 MALE_16 MALE_17 MALE_18 MALE_19 MALE_20_TO MALE_25_TO MALE_30_TO MALE_35_TO MALE_40_TO MALE_45_TO MALE_50_TO MALE_55_TO MALE_60_TO MALE_65_TO MALE_70_TO MALE_75_TO MALE_80_TO MALE_85 MALE_MEDIA FEMALE_ALL FEMALE_4 FEMALE_5_T FEMALE_10_ FEMALE_15_ FEMALE_15 FEMALE_16 FEMALE_17 FEMALE_18 FEMALE_19 FEMALE_20_ FEMALE_25_ FEMALE_30_ FEMALE_35_ FEMALE_40_ FEMALE_45_ FEMALE_50_ FEMALE_55_ FEMALE_60_ FEMALE_65_ FEMALE_70_ FEMALE_75_ FEMALE_80_ FEMALE_85 FEMALE_MED MARRIED_AG MARRIED_OR MARRIED COMMON_LAW UNMARRIED SINGLE SEPARATED DIVORCED WIDOWED MARRIED_A1 MARRIED_O1 MARRIED_M COMMON_LA1 UNMARRIED_ SINGLE_M SEPARATED_ DIVORCED_M WIDOWED_M MARRIED_A2 MARRIED_O2 MARRIED_F COMMON_LA2 UNMARRIED1 SINGLE_F SEPARATED1 DIVORCED_F WIDOWED_F FAMILIES_I FAMILY_SIZ FAMILY_SI1 FAMILY_SI2 FAMILY_SI3 COUPLE_FAM COUPLE_MAR COUPLE_MA1 COUPLE_MA2 COUPLE_MA3 COUPLE_MA4 COUPLE_MA5 COUPLE_COM COUPLE_CO1 COUPLE_CO2 COUPLE_CO3 COUPLE_CO4 COUPLE_CO5 SINGLE_PAR SINGLE_PA1 SINGLE_PA2 SINGLE_PA3 SINGLE_PA4 SINGLE_PA5 SINGLE_PA6 SINGLE_PA7 SINGLE_PA8 CHILDREN_F CHILDREN_1 CHILDREN_2 CHILDREN_3 CHILDREN_4 CHILDREN_5 POPULATIO1 POPULATIO2 POPULATIO3 POPULATIO4 POPULATIO5 POPULATIO6 POPULATIO7 POPULATIO8 POPULATIO9 POPULATI10 POPULATI11 POPULATI12 PRIVATE_HO PRIVATE_HH PRIVATE_H1 PRIVATE_H2 PRIVATE_H3 PRIVATE_H4 PRIVATE_H5 PRIVATE_H6 PRIVATE_H7 PRIVATE_H8 PRIVATE_H9 PRIVATE_10 PRIVATE_11 PRIVATE_12 PRIVATE_13 PRIVATE_14 PRIVATE_15 OCC_PRIVAT OCC_PRIVA1 OCC_PRIVA2 OCC_PRIVA3 OCC_PRIVA4 OCC_PRIVA5 OCC_PRIVA6 OCC_PRIVA7 OCC_PRIVA8 OCC_PRIVA9 PRIVATE_16 PRIVATE_17 PRIVATE_18 PRIVATE_19 PRIVATE_20 PRIVATE_21 PRIVATE_22 PRIVATE_23 NATIVE_LAN NATIVE_LA1 NATIVE_LA2 NATIVE_LA3 NATIVE_LA4 NATIVE_LA5 NATIVE_LA6 NATIVE_LA7 NATIVE_LA8 NATIVE_LA9 NATIVE_L10 NATIVE_L11 NATIVE_L12 NATIVE_L13 NATIVE_L14 NATIVE_L15 NATIVE_L16 NATIVE_L17 NATIVE_L18 NATIVE_L19 NATIVE_L20 NATIVE_L21 NATIVE_L22 NATIVE_L23 NATIVE_L24 NATIVE_L25 NATIVE_L26 NATIVE_L27 NATIVE_L28 NATIVE_L29 NATIVE_L30 NATIVE_L31 NATIVE_L32 NATIVE_L33 NATIVE_L34 NATIVE_L35 NATIVE_L36 NATIVE_L37 NATIVE_L38 NATIVE_L39 NATIVE_L40 NATIVE_L41 NATIVE_L42 NATIVE_L43 NATIVE_L44 NATIVE_L45 NATIVE_L46 NATIVE_L47 NATIVE_L48 NATIVE_L49 NATIVE_L50 NATIVE_L51 NATIVE_L52 NATIVE_L53 NATIVE_L54 NATIVE_L55 NATIVE_L56 NATIVE_L57 NATIVE_L58 NATIVE_L59 

You can obtain new (calculated) variables as follows. For instance, to obtain the proportion of residents who are between 20 and 34 years old, and between 35 and 49:

Hamilton_CT@data <- mutate(Hamilton_CT@data, Prop20to34 = (AGE_20_TO_ + AGE_25_TO_ + AGE_30_TO_)/POPULATION, Prop35to49 = (AGE_35_TO_ + AGE_40_TO_ + AGE_45_TO_)/POPULATION)

Example: Proportion of residents who are between 20 and 34 years old, and between 35 and 49:

Hamilton_CT@data <- dplyr::transmute(Hamilton_CT@data, AREA = AREA, TRACT = TRACT,
                                     POPULATION = POPULATION,
                                     POP20to34 = (AGE_20_TO_ + AGE_25_TO_ + AGE_30_TO_),
                                     Prop20to34 = POP20to34/POPULATION, 
                                     POP35to49 = (AGE_35_TO_ + AGE_40_TO_ + AGE_45_TO_),
                                     Prop35to49 = POP35to49/POPULATION, 
                                     POP50to64 = (AGE_50_TO_ + AGE_55_TO_ + AGE_60_TO_),
                                     Prop50to64 = POP50to64/POPULATION,
                                     POP65Plus = (AGE_65_TO_ + AGE_70_TO_ + AGE_75_TO_ + AGE_80_TO_ + AGE_85),
                                     Prop65Plus = POP65Plus/POPULATION)

This is a SpatialPolygonDataFrame. Convert to a dataframe (“tidy” it) for plotting using ggplot2:

Hamilton_CT.t <- tidy(Hamilton_CT, region = "TRACT")
Hamilton_CT.t <- rename(Hamilton_CT.t, TRACT = id)

Rejoin the data:

Hamilton_CT.t <- left_join(Hamilton_CT.t, Hamilton_CT@data, by = "TRACT")

Column `TRACT` joining character vector and factor, coercing into character vector

This is the function to create local Moran maps:

localmoran.map <- function(spat_pol = spat_pol, listw = listw, VAR = VAR, ID = ID){
  require(tidyverse)
  require(broom)
  require(spdep)
  require(plotly)
  
  spat_pol@data <- data.frame(ID = ID, VAR = VAR)
  spat_pol.t <- broom::tidy(spat_pol, region = "ID")
  spat_pol.t <- dplyr::rename(spat_pol.t, ID = id)
  spat_pol.t <- dplyr::left_join(spat_pol.t, spat_pol@data, by = "ID")
  
  df_msc <- transmute(spat_pol@data, 
                      ID = ID,
                      Z = (VAR-mean(VAR)) / var(VAR), 
                      SMA = lag.listw(listw, Z),
                      Type = factor(ifelse(Z < 0 & SMA < 0, "LL", 
                                           ifelse(Z > 0 & SMA > 0, "HH", "HL/LH"))))
  
  local_I <- localmoran(spat_pol$VAR, listw)
  
  spat_pol.t <- left_join(spat_pol.t, 
                             data.frame(ID = spat_pol$ID, local_I))
  spat_pol.t <- rename(spat_pol.t, p.val = Pr.z...0.)
  spat_pol.t <- left_join(spat_pol.t, 
                             df_msc)
  
  map <- ggplot(data = spat_pol.t, 
                aes(x = long, y = lat, group = group, 
                    p.val = p.val, VAR = VAR)) +
    geom_polygon(aes(fill = Type, color = p.val < 0.05)) +
    scale_fill_brewer(palette = "RdBu") +
    scale_color_manual(values = c(NA, "Black") ) +
    labs(color = "Prob < 0.05") +
    coord_equal() +
    theme(legend.title = element_blank())
  ggplotly(map, tooltip = c("p.val", "VAR"))
}

This is function is used to create \(G_i^*\) maps:

gistar.map <- function(spat_pol = spat_pol, listw = listw, VAR = VAR, ID = ID){
  require(tidyverse)
  require(broom)
  require(spdep)
  require(plotly)
  
  spat_pol@data <- data.frame(ID = ID, VAR = VAR)
  spat_pol.t <- broom::tidy(spat_pol, region = "ID")
  spat_pol.t <- dplyr::rename(spat_pol.t, ID = id)
  spat_pol.t <- dplyr::left_join(spat_pol.t, spat_pol@data, by = "ID")
  
  df.lg <- localG(VAR, listw)
  df.lg <- as.numeric(df.lg)
  df.lg <- data.frame(Gstar = df.lg, p.val = 2 * pnorm(abs(df.lg), lower.tail = FALSE))
  
  df.lg <- mutate(df.lg, 
              Type = factor(ifelse(Gstar < 0 & p.val <= 0.05, "Low Concentration",
                                   ifelse(Gstar > 0 & p.val <= 0.05, "High Concentration", "Not Signicant"))))
  spat_pol.t <- left_join(spat_pol.t,
                             data.frame(ID = spat_pol$ID, df.lg))
  map <- ggplot(data = spat_pol.t, 
                aes(x = long, y = lat, group = group, 
                    p.val = p.val, VAR = VAR)) +
    geom_polygon(aes(fill = Type, color = p.val < 0.05)) +
    scale_fill_brewer(palette = "RdBu") +
    scale_color_manual(values = c(NA, "Black") ) +
    labs(color = "Prob < 0.05") +
    coord_equal() +
    theme(legend.title = element_blank())
  ggplotly(map, tooltip = c("p.val", "VAR"))
}

Create spatial weights.

1. By contiguity:

Hamilton_CT.w <- nb2listw(poly2nb(pl = Hamilton_CT))

2. Binary, by distance (3 km threshold) including self.

Hamilton_CT.3knb <- Hamilton_CT %>% coordinates() %>% dnearneigh(d1 = 0, d2 = 3, longlat = TRUE)
Hamilton_CT.3kw <- nb2listw(include.self(Hamilton_CT.3knb), style = "B")

Activity

1. Create local Moran maps for the population and proportion of population in the age group 20-34. What is the difference between using population (absolute) and proportion of population (rate)? Is there a reason to prefer either variable in analysis? Discuss.

localmoran.map(Hamilton_CT, Hamilton_CT.w, Hamilton_CT$POP20to34, Hamilton_CT$TRACT)

Loading required package: plotly

Attaching package: <U+393C><U+3E31>plotly<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:ggplot2<U+393C><U+3E32>:

    last_plot

The following object is masked from <U+393C><U+3E31>package:stats<U+393C><U+3E32>:

    filter

The following object is masked from <U+393C><U+3E31>package:graphics<U+393C><U+3E32>:

    layout

Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`

localmoran.map(Hamilton_CT, Hamilton_CT.w, Hamilton_CT$Prop20to34, Hamilton_CT$TRACT)

Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`

2. Use the \(G_i^*\) statitic to analyze the population and proportion of population in the age group 20-34. What is the difference between using population (absolute) and proportion of population (rate)? Is there a reason to prefer either variable in analysis? Discuss.

gistar.map(Hamilton_CT, Hamilton_CT.3kw, Hamilton_CT$POP20to34, Hamilton_CT$TRACT)

Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`

gistar.map(Hamilton_CT, Hamilton_CT.3kw, Hamilton_CT$Prop20to34, Hamilton_CT$TRACT)

Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`

3. Now create local Moran maps for the population and population density in the age group 20-34. What is the difference between using population (absolute) and population density (rate)?

localmoran.map(Hamilton_CT, Hamilton_CT.w, Hamilton_CT$POP20to34, Hamilton_CT$TRACT)

Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`

localmoran.map(Hamilton_CT, Hamilton_CT.w, Hamilton_CT$POP20to34/Hamilton_CT$AREA, Hamilton_CT$TRACT)

Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorJoining, by = "ID"
Column `ID` joining character vector and factor, coercing into character vectorWe recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`

4. More generally, what do you think should guide the decision of whether to analyze variables as absolute values or rates?